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ABSTRACT 

We present WHFast, a fast and accurate implementation of a Wisdom-Holman symplectic 
integrator for long-term orbit integrations of planetary systems. WHFast is significantly faster 
and conserves energy better than all other Wisdom-Holman integrators tested. 

We achieve this by significantly improving the Kepler-solver and ensuring numerical sta¬ 
bility of coordinate transformations to and from Jacobi coordinates. These refinements allow 
us to remove the linear secular trend in the energy error that is present in other implementa¬ 
tions. For small enough timesteps we achieve Brouwer’s law, i.e. the energy error is dominated 
by an unbiased random walk due to fioating-point round-off errors. 

We implement symplectic correctors up to order eleven that significantly reduce the en¬ 
ergy error. We also implement a symplectic tangent map for the variational equations. This 
allows us to efficiently calculate two widely used chaos indicators the Lyapunov characteristic 
number (LCN) and the Mean Exponential Growth factor of Nearby Orbits (MEGNO). 

WHFast is freely available as a flexible C package, as a shared library, and as an easy-to- 
use python module. 

Key words: methods: numerical — gravitation — planets and satellites: dynamical evolution 
and stability 


1 INTRODUCTION 

Celestial mechanics, the field that deals with the motion of celes¬ 
tial objects, has been an active field of research since the days of 
Newton and Kepler. Analytic solutions only exist for a few special 
cases. Historically, the main driver for the development of pertur¬ 
bation theory has been the problem of planets orbiting the Sun. Be¬ 
cause the central body is so much more massive than the planets, it 
is profitable to ask how the small mutual tugs between the planets 
modify the Keplerian orbits they would each individually follow 
around the Sun in the absence of the other bodies. This analytical 
approach has been, and continues to be, successful in explaining 
many important features of planetary orbits. However, the Solar 
System is chaotic, and the rise of computing power has yielded 
many important insights. There is therefore considerable interest in 
developing fast and accurate numerical integrators. 

A large number of such integrators have been developed over 
the years to perform this task. For many long term integrations, 
symplectic integrators have proven to be a favourable choice. Sym¬ 
plectic schemes incorporate the symmetries of Hamiltonian sys¬ 
tems, and therefore typically conserve quantities like the energy 
and angular momentum better than non-symplectic integrators. 


For integrations of planetary systems, [Wisdom & Holm^ 
( |1991| ), and independently [Kinoshita et ah] ( 199l| ), developed a 
widely used class of symplectic integrators. The ideas of jWisdom] 
|& Holman|(|1991jl de veloped from the original ideas of the mapping 
method of |Wisdom] ( |1981| ). We refer to these as a Wisdom-Holman 
mapping or a Wisdom-Holman integrator. Since then, many authors 
have modified and built upon this method, and several have made 
their integrators publicly available to the astrophysics community 
(e.g. [Chambers & Migliorinijl 997 [ [Duncan et al.|1998| l. 


The Wisdom-Holman integrator exploits the intuition from 
perturbation theory that one can separate the problem into a system 
of Keplerian orbits about the Sun, modified by small perturbations 
among the planets. The nuisance is that while Newton provided 
us the solution to the two-body problem, Poincare showed that the 
remaining superimposed perturbations are not integrable. Analyt¬ 
ically, the traditional way forward is to average over the short- 
period oscillations in the problem to yield approximate solutions. 
The great insight of Wisdom and Holman was that, at the same level 
of approximation, one can add high frequency terms. By judicious 
choice of these additional frequencies, the perturbations among the 
planets can be transformed into trivially integrated delta functions. 
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The result is an exceedingly efficient integrator that has proven an 
indispensable tool for modern studies in celestial mechanics. 

In this paper, we present results from a complete reimple¬ 
mentation of the Wisdom-Holman integrator. We show how to 
speed up the algorithm in several ways and dramatically increase 
its accuracy. Many of the improvements are related to finite dou¬ 
ble floating-point precision on modern computers (IEEE754, |ISO| 
|2011| l. The fact that almost all real numbers cannot be represented 
exactly in floating-point precision leads to important consequences 
for the numerical stability of any algorithm and the growth of nu¬ 
merical round-off error. 

To our knowledge, we present the first publicly available 
Wisdom-Holman integrator that is unbiased, i.e. the errors are ran¬ 
dom and uncorrelated. This leads to a very slow error growth. Eor 
sufficiently small timesteps, we achieve Brouwer’s law, i.e., the en¬ 
ergy error grows as time to the power of one half. 

We have also sped up the integrator through various improve¬ 
ments to the integrator’s Kepler-solver. Our implementation allows 
for the evolution of variational equations (to determine whether or¬ 
bits are chaotic) at almost no additional cost. Additionally, we im¬ 
plement so-called symplectic correctors up to order eleven to in¬ 
crease the accurary ( [Wisdom et al.|1996] l, allow for arbitrary unit 
choices, and do not tie the integration to a particular frame of ref¬ 
erence. 

We make our integrator, which we call WHFast, publicly avail¬ 
able in its native C99 implementation and as an easy-to-use python 
module. 

The remainder of this paper is structured as follows. We first 
summarize the concepts and algorithms used in this paper, includ¬ 
ing Jacobi Coordinates, our choice of Hamiltonian splitting, the 
symplectic Wisdom-Holman map, symplectic correctors and the 
variational equations in Sect. We then go into detail discussing 
the improvements we have made to these algorithms in Sect.j^ Nu¬ 
merical tests are presented in Sect. [^before we conclude in Sect.|^ 


2 BACKGROUND 

The Hamiltonian TT of the gravitational A-body system can be writ¬ 
ten as the sum of kinetic and potential terms in Cartesian coordi- 


nates 
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One way forward toward separating out the two-body Keplerian 
Hamiltonians is to transform to heliocentric coordinates involving 
the centre-of-mass and the r/ -Tq. However, rewriting the Cartesian 
momenta in terms of heliocentric momenta (which have an addi¬ 
tional component along the centre-of-mass momentum), leads to 
several cross-terms. Alternatively, Jacobi worked out a coordinate 
system in which the kinetic terms are particularly clean, and the ki¬ 
netic energy remains a sum of squares. Eor readers that may not be 
familiar, and because our improved accuracy is largely due to mod¬ 
ifications of the manner in which we transform between Cartesian 
and Jacobi coordinates, we briefl y review them (see alsojPlummerj 
|1918[[Sirssman & Wisdom|2001 [[Murray & Dermott|2000|~ 


2.1 Jacobi Coordinates 


Rather than reference planet positions to the central star, a planet’s 
Jacobi coordinates are measured relative to the centre-of-mass of 
all bodies with lower indices. Eor concreteness, consider a system 
of N particles with masses m/, / = 0,..., A - 1. Let r/ be the posi¬ 
tion vector of the i-th particle with respect to an arbitrary origin that 
is fixed in an inertial frame. Here we assume that the particles are 
ordered such that i = 0 corresponds to the central object, / = 1 to 
the innermost object orbiting the central object and so on. The ex¬ 
istence of such an ordering does not restrict the architecture of the 
system. Eor example, the coordinates of an equal-mass binary with 
a circumbinary particle can be expressed in Jacobi coordinates. But 
note that the ordering might in general be non-unique and that it can 
change during an integration. This can have important implications 
for a numerical scheme using Jacobi coordinates. 

The Jacobi coordinate rj of the /-th particles is the position 
relative to R/_i, the centre-of-mass of all the particles interior to 
the i-th particle: 

r- = Yi - R/_i, for i = 

1 * 

where R, = — > m/r, and M, = 

Mi ^ ^ 


1,...,A-1 (2) 

i 

J]mj. (3) 
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Other quantities such as the velocity and acceleration (also the co¬ 
ordinates in the variational equations, see below) transform in the 
same way. This is because the Jacobi coordinates are a linear func¬ 
tion of the Cartesian coordinates, and the velocity is the time deriva¬ 
tive of the position in both coordinates systems. 

The momenta, however, transform differentl}Q The momen¬ 
tum conjugate to r^' and the corresponding Jacobi mass are given 
by 


pj = m-r'i = rn'iVi 


and 


Mi 


nii Mi_i 
rui + M/_i ’ 


( 4 ) 


Note that the Jacobi mass m'. is the reduced mass of m/ and M/_i. 
Explicit expressions for the momenta can be found by evaluating 
the time derivative rEq.|^ 

The Jacobi coordinates above are relative coordinates for i = 
1,..., A - 1. Eor the 0-th coordinate, a different convention is used. 


N-l 

To = Ra,_i, m'a = Mn-i, Po = 

>0 

Thus, Fq points towards the centre-of-mass of the entire sys¬ 
tem, Pq is the total momentum and is the total mass. 


2.2 Hamiltonian Splitting 

After some algebra, we can rewrite the Hamiltonian in Eq. in 
terms of the conjugate momenta of the Jacobi coordinates (e.g. 
[Murray & Dermott|[2000[ [Sussman & Wisdom[[2001[ l. We only 
rewrite the kinetic term and keep the potential term expressed as 
a function of the Cartesian coordinates: 
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^ But note that we do not need to calculate the momenta explicitly in our 
algorithm. 
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Note that the kinetic term is still diagonal, i.e. there are no cross 
terms involving p/py with i ^ j. Next, we add and subtract the term 




Z 


Gm\Mi 

|r;i ■ 


(7) 


After grouping terms in the Hamiltonian, we arrive at 
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*^Interaction 

The first term, *770, simply describes the motion of the centre-of- 
mass Tq along a straight line. For that reason this term is often ig¬ 
nored. However, we keep it which will allow us to integrate parti¬ 
cles without any restriction to a particular frame of reference. 

The terms TTKepier can be split up further into a sum of 


^TTKepler^^. 


pf Gm'M; 
2m: |r;| 


(9) 


Each of the Hamiltonians ('^Kepler), describes the Keplerian motion 
of the /-th particle with mass m/ around the centre-of-mass of all 
interior particles with total mass M/_i. 

After some more algebra, the interaction term can be simpli¬ 
fied and split into two parts, one of which can be easily computed 
in Jacobi coordinates and the other in Cartesian coordinates of the 
inertial frame 
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One important point to note is that our choice of *77+ is slightly 
different from that used by [Murray & Dermott| ( |2000j ) and |Wisdom| 
|& Holman| ( |1991| ). These authors use 


( 77 +) wh 1991 


z 


Gm'. Mi 


(11) 


where Mi = Their choice leads to the usual disturbing 

function in perturbation theory. We conducted various tests but 
found no significant difference between these mass choices. We 
therefore chose our prescription, Eq.|^ which has a simpler phys¬ 
ical interpretation: the mass entering Kepler’s third law is simply 
the interior mass. 


2.3 Wisdom-Holman Mapping 

Our goal is to find a solution to the equations of motion for par¬ 
ticles governed by the Hamiltonian in Eq. No analytic solution 
exists to the full Hamiltonian and we thus need to find an approxi¬ 
mate solution. There are many different ways to do that. Here, we 
describe the idea of constructing a symplectic integrator by means 
of splitting the Hamiltonian into smaller parts, each of which can 
be easily integrated. 

The introduction of Jacobi coordinates led us to the Hamil¬ 
tonian splitting described in Sect. |2.2| Analytic solutions can be 
found for the evolution of the system under each of the individual 


Hamiltonians TTo and TTinteraction- The solution to TTo simply corre¬ 
sponds to motion along a straight line. The solution to ‘77interaction is 
a kick step where the velocities change due the inter-particle accel¬ 
erations but the position remain constant. The solution to TTKepier 
is a set of two-body Kepler orbits, which can also be easily solved 
with an iterative algorithm. We discuss the details related to the 
Kepler problem in Sect. |2.6| 

Now that we have broken down the full Hamiltonian into in¬ 
dividual Hamiltonians, to all of which we know the solution (or 
can easily calculate it), we can construct a symplectic integrator for 
the total Hamiltonian using an operator split method (e.g. [Saha &| 
|Tremaine]|1992| ). Let us describe the evolution of particles under 
a Hamiltonian *77 for a time dt using the operator notation 
The notation ^ 2 {dt) o i^i(dt) means applying operator ‘Ki first, 
then applying operator ‘Hi- It is easy to see that many of the oper¬ 
ators commute with each other, i.e. 


[tVo, '^Kepler] - 0 

(12) 

[TTq, *77interactionj “ 0 

(13) 

"H. 

II 

O 

< 

(14) 

where [TVi , ^ 2 ] = “Hi o tVz - tVz o -Hi 

. This leads to the following 


Drift-Kick-Drift (DKD) operator splitting scheme, which we refer 
to as the Wisdom-Holman map: 

(Drift) Evolve the system under ‘KKepier(<7t/2) o i^Q(dt/2). 

(Kick) Evolve the system under ‘Hinteraction(<70- 

(Drift) Evolve the system under i^Kopicridt/2) o i^Q(dt/2). 

The ordering of ‘KKepier and TTo in the first and last step doesn’t 
matter as they commute. The first and last steps can be combined if 
the system is evolved for multiple timesteps. 

Note that the evolution of TTKepier and TTo is most easily 
accomplished in Jacobi coordinates. The interaction Hamiltonian 
77interaction, howcvcr. Contains terms that depend on both the Carte¬ 
sian and Jacobi coordinates. The simplest way to calculate these 
terms is to convert to Cartesian coordinates, evaluate the r/ - Vj 
term, convert the accelerations back to Jacobi accelerations, and 
calculate the remaining terms. 


2.4 Symplectic Correctors 


The operator splitting method used in the symplectic integrator dis¬ 
cussed above effectively adds high frequency terms to the Hamil¬ 
tonian. An argument often used in favour of symplectic integrators 
is that, although these high-frequency terms alter the Hamiltonian, 
they do not change the long term evolution as they average out. 
However, they do lead to relatively large short term oscillations, 
for example in the energy error. 

The idea of a symplectic corrector, first used by Tittemore"^ 
|Wisdom| ( |1989| ) and fully developed by [Wisdom et al.| ( |1996| l, is 
to remove some of these high frequency terms using perturbation 
theory. The basic procedure is as follows. Before the start of an in¬ 
tegration, we convert from real coordinates to so-called mapping 
coordinates. Then we perform the integration using our standard 
symplectic map. After the simulation has finished (or whenever we 
need an output) we convert back from mapping to real coordinates. 
The symplectic corrector operator that we use is a combination of 
several ‘Hinteraction(<70 and i^KtpiAdt) o iHo{dt) operators applied 
for different (positive and negative) intervals dt. If e is the order 
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of the perturbations, i.e. the mass ratio and therefore the relative 
magnitude of ‘Kinteraction compared to ‘KKepier, then one can show 
that the use of symplectic correctors can lead to a scheme of order 
0{edt^) + 0(6^df) where K is the order of the symplectic corrector 
( [Mikkola & Palmer|2000| l. A second order Wisdom-Holman map 
without symplectic correctors has an energy error of order 0{e^dt^). 
Because this coordinate transformation for the symplectic corrector 
is only performed for outputs and at the beginning and end of the 
simulation, its effect on the speed of the algorithm is negligible for 
sparse output. 

A full derivation of the symplectic correctors would go beyond 
the scope of this paper and we refer the reader to [Wisdom et al.| 
\\996) and [Mikkola & Palmer] ( [2000[ l. The corrector coefficients 
are listed in a compact form in Wisdom] { |2006[ ). 

We implement a third, fifth, seventh and eleventh order sym¬ 
plectic corrector for WHFast. Whether the high order-symplectic 
correctors provide any improvement over the low-order ones de¬ 
pends on the mass ratios in the system. For Jupiter-mass planets, a 
symplectic corrector of fifth order is no less accurate than a higher 
order one. If in doubt, there is no harm done in using a higher- 
order corrector as the speed implications are minimal. Thus, we 
implement the eleventh-order symplectic corrector by default. 


2.5 Chaos Indicators 

A powerful tool for studying the long term evolution of Hamil¬ 
tonian systems is the Lyapunov characteristic number (LCN). The 
inverse of the LCN is the Lyapunov timescale and gives an estimate 
of how fast two nearby particle trajectories diverge. If the system 
is chaotic, the divergence is exponential in time and the Lyapunov 
timescale is finite. Thus, measuring the LCN gives us an estimate 
of whether the system is chaotic and, if so, on what timescale. 

A more recent approach with similar informative value is the 
Mean Exponential Growth factor of Nearby Orbits, or MEGNO for 
short ( [Cincotta et al.[2003] l. The MEGNO, Y{t), is a scalar function 
of time, and provides a clear picture of resonant structures and of 
the locations of stable and unstable periodic orbits. 

There are two ways to calculate the LCN or the MEGNO. 
Conceptually the simplest is to integrate an additional shadow par¬ 
ticle for each body in the simulation, i.e. a particle with slightly per¬ 
turbed initial conditions. One can then directly measure the diver¬ 
gence of each particle’s path from its shadow. The second approach 
is to consider each body’s six-dimensional displacement vector 6i 
from its shadow (in both position and velocity) as a dynamical vari¬ 
able. One can then obtain differential equations for each 6i vector 
by applying a variational principle to the trajectories of the origi¬ 
nal bodies. We choose to follow the latter approach, as it is both 
faster and numerically more robust ( [Tancredi et al.|20()T] ). In this 
scheme, one can imagine shadow particles with phase-space coor¬ 
dinates +Si, where= (ri, Vi) is the phase-space coordinate 

of the i-th original particle. Initially, we set each component of 6i 
to a small value. 

We follow the work of [Mikkola & Innanen] ( [19991 ) who de¬ 
scribe how to efficiently couple the variational equations to the 
original equations of motion. This allows us to construct a sym¬ 
plectic integrator for the variational equations (a symplectic tangent 
map). An important advantage of this method is that we only solve 
Kepler’s equation once for each particle/shadow-particle pair (one 
of the most time-consuming steps in a Wisdom-Holman integrator 
for small particle numbers). 


The MEGNO is then straightforwardly computed from the 
variations as ( [Cincotta et al.[2003] l 


Y(t) = 


2 

tJo ZZo^f(t') 


(15) 


If Y(t) oo, then the system is chaotic. For quasi-periodic orbits, 
the MEGNO converges to a finite value, Y(t) 2 (e.g.[Hinse et al.[ 

[20T0l ). 

One can obtain the Lyapunov characteristic number (LCN), 
the inverse of the Lyapunov timescale, from the time evolution of 
the MEGNO via a linear least square fit to Y(t). 


2.6 Kepler Problem with Variations 

In this section, we summarize how to solve the two-body Kepler 
problem numerically, including the variational equations. Although 
the solution has been known since the days of Newton, the tran¬ 
scendental nature of Kepler’s equation does not admit a closed- 
form mathematical expression. 

We closely follow the work of ( [Mikkola & Innanen[[1999] ) 
where the reader can find additional information that we have left 
out. The equivalent one-body Hamiltonian for the Kepler problem 
is 

^^Kepler = 

where M is the total mass of the two bodies. For consistency with 
[Mikkola & Innan^ \\999) , we have dropped the primes, have 
scaled out from p'., and rewritten Eq.j^in non-dimensional form, 
i.e. the gravitational constant G = 1 for the remainder of this pa¬ 
per. However, we have taken care to remove any dependence on the 
choice of units from our implementation, so G can be freely set by 
the user in our implementation of the algorithm. 

Our task is to find the final positions and velocities r and v of a 
particle evolving under this Hamiltonian for some time dt, given the 
initial conditions Vq and Vq. Thus, we seek the effect of the operator 

^Kepler (^0- 

It is advantageous to solve the Kepler problem numerically 
using the Gauss f and g functions, which express the relevant quan¬ 
tities in terms of Tq and Vq ( [Wisdom & Holman] 199 1[ ). This avoids 
the computationally expensive conversion between Cartesian and 
classical orbital elements, and avoids coordinate singularities asso¬ 
ciated with circular orbits. We find that it is advantageous to use 
universal variables in this solution ( [Stumpff|1962[ ). This approach 
provides greater speed and numerical stability compared to a solu¬ 
tion using elliptic elements. It also avoids the singularity associated 
with the transition from elliptic to hyperbolic motion. 

To solve the analogue of Kepler’s equation for the particle’s 
position in time, we make use of several special functions. Let us 
begin by defining the c-functions ( [Stumpff[1962[ ) as a series expan¬ 
sion: 


oo 


Cn(z) = L 

y=0 


(-zy 
(n + 2jy: 


which satisfy the recursion relation 

Cni.^ ~ T ~ ZCn+2’ 
n\ 


(17) 


(18) 
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The c-functions are related to trigonometric functions, for example 

sin yi 

co(z) = cos ^Jz and ci(z) = -(19) 

^/z 

and thus satisfy the following relationships ( |Mikkola|1997| ), which 
are related to the half-angle formula for trigonometric functions: 

C5(z) = ^[C5(z/4) + C4(z/4)+C3(z/4)C2(zm (20) 

C4(z) = ^C3 (z/4)[1+ci(z/4).] (21) 

Values for cq through C 3 are then readily co mputed from Eg. 
Next, we introduce the so called G-functions ( [Stiefel & Scheifelq 
|1971| l which in turn depend on the c-functions: 

Gn{p.X) = X^Cn{px\ (22) 

The G-functions also satisfy recursion relationships similar to 

those mentioned above for the c-functions. We can easily calcu¬ 
late derivatives of G„ by looking at the series expansion of c„ (see 
[Mikkola & Innanen|1999| for details). With this framework, we can 
now write down the steps needed to find the solution to the Kepler 
Hamiltonian in compact form. 

First, we need to calculate the following three quantities from 


the 

initial conditions Tq, Vq: 



2M , 


P 

- Vq 

(23) 


ro 


m 

= ro • Vo 

(24) 

(0 

= M-pro 

(25) 


where ro = |ro| and vq = |vo|. Note that the semi-major axis a can 
be written as a = M//5. 

Second, we need to solve Kepler’s equation which, using the 
above notation, takes the form 


roX -h rioGiiP. X) + X) - dt ^ 0. 


(26) 


We solve this equation for X. This is a non-algebraic (i.e. transcen¬ 
dental) equation that we need to solve iteratively, for example using 
Newton’s method. In Sect. |3.2| we describe our algorithm in detail. 

Third, having solved Kepler’s Equation, we can calculate the 
so called GauB / and g-functions as well as their time derivatives 
via 


f^l-M 


G 2 

ro 


g^dt- MG 3 


MGi 

ror 

MG2 


(27) 

(28) 


where r = ro rjoGi + ^ 0 ^ 2 - Note that all the G-functions depend 
on p and the X value found in the second step. 

Fourth, we write the final positions and velocities as a linear 
transformation of the initial conditions using the GauB / and g- 
functions: 


r = /ro + gvo V = /ro + gvo. (29) 

This completes the solution of the Kepler problem. 

To solve for the variational equations, we also make use of the 
G-functions. Fortunately, we only need to solve Kepler’s equation 
once (to solve for X). We then get the solution for the variational 
equations without solving another transcendental equation and thus 
have only one iteration loop per timestep for both the particle and 
its variational counterpart. The position and velocity components 


of 6 at the end of the timestep, and ^v, can be written as 


6 r = f 6 ro + g 6 \o + ro 6 f + Yo 6 g (30) 

= /^ro+ g^Vo+ ro ^/ +Vo ^g, (31) 


where the variations Sf, Sg, 6 f and Sg can be derived from Eqs. 27 
(see [Mikkola & Innanen|1999| for the explicit expressions). 


2.7 Types of Numerical Errors 

There are three distinct effects contributing to the energy error of a 
symplectic integrator (see e.g., [Quinn & Tremaine|199()] l. See also 
[Rein & Spiegel| ( [^15| ) for a similar discussion for non-symplectic 
integrators. 

First, there is an error term associated with the integrator it¬ 
self because we are not solving the equations of motion for the 
Hamiltonian 7/ exactly. For symplectic integrators such as those 
discussed here, this error term is bound and we call it Fbound- If the 
mass ratio of the planets to the star is 6 , then the order of this error 
term is roughly 0 {edt^) for integrators without symplectic corec¬ 
tors and 0(6 dt^) + 0 (e^ dt^) for those with symplectic correctors 
(see Sect. [2.4j ). Note that Ebound is independent of time t. 

Second, there is an error term associated with the finite preci¬ 
sion of numbers represented on a computer. We can only represent a 
small subset of all real numbers exactly in floating-point precision. 
Thus after every operation such as an addition or multiplication, 
the computer rounds to a nearby floating-point number. For CPUs 
and compilers that follow the IEEE754 standard ( [ISO|201l'[ ), we are 
guaranteed to round to the nearest floating-point number. Thus, if 
all operations follow the IEEE754 standard, then as long as the al¬ 
gorithm itself is unbiased, we expect the error to grow as the square 
root of the number of operations, i.e. Erand ~ ~ where N 

is the number of timesteps. This is the best behaviour achievable; 
to do better we would have to move to extended precision or use 
fewer operations. This fundamental limit is known as Brouwer’s 
law ( [Newcombj 1899 [[Brouwer 1 19T7] l. 

Third, if any parts of the integration algorithm are biased, the 
errors will be correlated. This leads to a faster long-term energy- 
error growth than if errors are uncorrelated; it grows linearly with 
time, i.e. Fbias ~ ~ 

For a given integrator, which of these three error terms domi¬ 
nates depends on the nature of the simulation, the timestep, and the 
total integration time (number of timesteps). 


3 IMPROVEMENTS 

The algorithms we describe in Sect.[^have been used successfully 
for many years. In the following, we show how to significantly im¬ 
prove the speed and accuracy of the algorithms by taking special 
care in the implementation of several details, many of which are 
related to finite floating-point precision on modern computers. 

For the remainder of this paper, we will assume that we work 
with a CPU that follows the IEEE 754 standard for floating-point 
arithmetic. Most importantly, we assume that all floating-point op¬ 
erations follow the rounding to nearest, ties to even rule ( [ISO[ 
[2011| ). What follows is in principle applicable to any precision. 
However, we work exclusively in double floating-point precision 
(64 bit) which is used on almost all modern CPUs. 


© 0000 RAS, MNRAS 000, 000-000 


























6 


Hanno Rein and Daniel Tamayo 


3.1 Jacobi Coordinate Transformations 

The evolution under the effect of the interaction Hamiltonian is 
most efficiently done in Cartesian coordinates. On the other hand, 
the evolution of the Kepler Hamiltonian is easier in Jacobi coordi¬ 
nates. We thus need an efficient way to convert to and from Jacobi 
coordinates. 

Luckily, the conversion from Cartesian to Jacobi coordinates 
and back can be done efficiently in 0{N). We construct the algo¬ 
rithms from the definitions above and list them here in pseudo code. 
As before, primes denote Jacobi coordinates. Note that these algo¬ 
rithms work even if some of the bodies are test particles with m/ = 0 
(for / ^ 0). To convert from Cartesian to Jacobi coordinates: 

R ^ mo • ro 
for / <— 1, A - 1 do 
r; ^ Yi - 

R <— R • (1 -h + rui • r'. 

Tq <— RIMn_i > This is the centre-of-mass. 

Similarly, we construct the algorithm to convert back from Jacobi 
to Cartesian coordinates as follows: 

R <— Fq • > Centre of mass, 

for / <— A - 1,1 do > Loop is in reverse order. 

R^(R-mrr;)/M, 

Yi ^ r; -h R 
R ^ R • Mi_i 

ro <— R/mo > Setting the coordinate of the 0-th particle. 

We thoroughly tested the conversions to and from Jacobi coordi¬ 
nates to ensure they are unbiased. This task turns out to be much 
harder than we naively expected. As an example, consider the fol¬ 
lowing algorithm which is formally equivalent to the above but nu¬ 
merically much less stable. 

R^O 

for / <— A - 1,1 do > Loop is in reverse order. 

Yi <— Fq -h • rj - R I’d is the centre-of-mass. 

R <— R -h m/ jMi • y\ 

Fq <— Fq - R > Setting the coordinate of the 0-th particle. 

In the above algorithm, we access r^ multiple times and have to 
do a subtraction in the last step. This significantly promotes error 
propagation and leads to floating-point errors that can be orders of 
magnitudes higher than in the other implementation. After many 
timesteps, this leads to a linear secular growth in the energy error. 


3.2 Implementation of Newton’s Method 

To solve Kepler’s equation (Eq.[^ for X, we need to use an iter¬ 
ative scheme. We now describe our implementation of Newton’s 
method in floating-point arithmetic. The straightforward imple¬ 
mentation is an iteration loop that terminates when the change to 
X is small, e.g., 

X <— initial guess 

Fepeat 

dX ^-f{X)/f{X) 

X^X + dX 
until \dX/X\ < eps. 


Here, eps is a small number just above machine precision, typi¬ 
cally eps ~ 10“^^. We use a different implementation of Newton’s 
method that is both faster and more accurate, despite the fact that it 
is algebraically equivalent to the above implementation. 

X <— initial guess 

Xprevi NaN > Any number different from X works. 

repeat 

-^prev2 ^ -^prevl 
^prevl ^ ^ 

X^{X-f{X)-f{X))/f{X) 
until X — Xpj-eyj or X — Xpj-ey 2 

Note that the equal sign in the above breakout condition is evaluated 
in floating-point precision. In comparison to the first algorithm, at 
each iteration step we test whether the iteration has converged by 
a simple comparison rather than by a slow division and absolute- 
value operation. 

We keep track of two previous values instead of just one be¬ 
cause for certain initial conditions, the iteration can cycle indefi¬ 
nitely between two nearby floating-point numbers and not converge 
to a single floating-point number. 

Our implementation thus ensures that the value of X is more 
accurately calculated than in the straightforward implementation 
using a heuristic value of eps. A further advantage of rewriting 
Newton’s method in the above form is that the term on the right- 
hand-side of the last line can be simplified significantly for the Ke¬ 
pler problem, giving: 

^ ^ + ^qG2) - T]oG2 - ^0^3 + dt 

ro rjoGi + (0G2 

where the G’s on the right-hand-side all depend on X and p (see 

Eq.|^. 

We also experimented with higher-order generalizations of 
Newton’s method (Householder’s methods). For typical cases 
where the orbits are not extremely elliptical {e < 0.99) and the 
timestep is much smaller than the shortest orbital period, we found 
Newton’s method to always be fastest. This is because when the 
value and derivatives of the function are easily evaluated, the pre¬ 
cision gain from these higher-order methods does not compensate 
for the increased computation cost of each iteration. In other words, 
while higher-order methods will converge in fewer iterations than 
Newton’s method, the overall computation time is longer. At large 
eccentricities and long timesteps, the G-function evaluations be¬ 
come expensive (one must recursively apply the quarter-angle for¬ 
mulas described in Sect. |3.5| ), and higher-order methods are help¬ 
ful. For large eccentricities we use a higher order method described 
in detail in Sect. |3.4| To safeguard against rare cases where New¬ 
ton’s method might fail, we also implemented a failsafe bisection 
method. We find that the bisection method is only triggered when 
the timestep is comparable to the orbital period. 

3.3 The Initial Guess foF KepleF’s Equation: ShoFt Timesteps 

The quantity X in Fq.j^can also be expressed as 

r^o+dt 7 f 

—=dt’(r-^) (33) 

Jto ^ 

where to is the time at the beginning of the timestep, and (r~^} is 
the time-averaged value of over the interval [to, to + dt]. Thus, 
if the orbit’s eccentricity e is low, or more generally if the timestep 
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is short enough that the orbital radius does not vary much, then 
X ^ dtlvQ. The troublesome cases are highly eccentric orbits near 
pericentre where the radius changes rapidly. For such cases, the 
radius varies by a factor of 1 + ^ 2 from pericentre to a true 

anomaly of 90°. We can therefore estimate the timescale over which 
the orbital radius varies near pericentre as 


Tchar — 


q 


a(\-e) l \-e '^'^ (1 - 

na \l + e) n 


(34) 


where q is the pericentre distance, is the speed at pericentre and 
n is the mean motion. Thus, if one does not resolve pericentre pas¬ 
sages (i.e., ndt = AM > (1 - X will differ from dtjrQ near 

pericentre (but may nevertheless conform to the simple approxima¬ 
tion at apocentre where the body moves slowly). 

More quantitatively, one can non-dimensionalize Eq. set¬ 
ting X - roXIdt. One can then solve the equation perturbatively, 
assuming the deviations from X - I are small. This procedure re¬ 
quires that the following three non-dimensional parameters in the 
equation also be much smaller than unity. 




~r^ 


qodt 


iodt^ 


(35) 


One can show that when our heuristic estimate AM (1 - is 
satisfied, ^ ^ 1 • In this case, one can extend the solution of 
Eq.j^to higher order. Eor the initial guess in our algorithm, we go 
up to second order 



(36) 


We experimented with higher-order initial guesses (see |Danby| 
|1987| for explicit expressions), but found these to be slower, even 
for small eccentricities and timesteps. This can again be attributed 
to the computational efficiency of each iteration of Newton’s 
method. 


3.4 Large Eccentricities and Timesteps 

The previous two sections describe an optimized algorithm for 
solving Kepler’s equation when the timestep and eccentricities 
are low. We have also developed an improved handling of high- 
eccentricity/long-timestep cases. In this regime, both the solver and 
initial guess should be modified. 

Like previous authors ( |Conway|1986||Danby|1987| ), we found 
the root-finding method of Laguerre-Conway to be most stable. 
However, unlike |Danby| ( |1987| ), who finds the method to always 
converge (presumably using comparatively small timesteps), we of¬ 
ten have to resort to bisection when the timestep is comparable to 
the orbital period. Of course, such long timesteps should not be 
chosen anyway, since they poorly sample inter-planet interactions, 
and are more susceptible to timeste p resonances ([Wisdom & Hol-| 
|man|1992l|Touma & Wisdom] 1993 [|^uch & Holman] 199^ 

We also had to modify the breakout condition used for New¬ 
ton’s method. While the Laguerre-Conway algorithm sometimes 
also bounces between two floating-point values once it has con¬ 
verged, in this regime the method often executes larger-period cy¬ 
cles (e.g., it will periodically repeat the last eight floating-point 
numbers). We therefore chose to store the values from each iter¬ 
ation and exit the loop whenever a result was repeated. 

One way to determine which solver should be used is to check 
whether dt is smaller than Tchar (Eq. [^. However, because Tchar 


is expensive to compute from Fq and Vq, we instead check how 
much the first iteration of Newton’s method deviates from the ini¬ 
tial guess, as a fraction of The latter is a natural quantity 

to compare against since it is the value of X when the timestep is 
equal to the orbital period. We found a threshold of ~ 1% to strike a 
good balance over a wide parameter range in timestep/eccentricity 
space, though the algorithm’s speed is not particularly sensitive to 
the exact value adopted. 

Einally, the method can be sped up in this regime with an 
improved initial guess for X, since Jt/ rp in Eq.[^blows u p near 
pericentre as the eccentricity gets large. [Danby & Burkardt| ( |1983| ) 
provide a widely used initial guess using classical orbital elements 
but, to our knowledge, no comparably simple initial guess has been 
found for universal variables. 

In this high-eccentricity / long timestep regime, most existing 
methods using universal variables choose to make the expensive 
conversion to orbital elements and use Danby’s guess. We instead 
observed that because (r"^) = a~^ over one orbital period, X = dt/a 
for a timestep of one orbit. We find that over a relevant parameter 
range with timesteps logarithmically spaced between 0.03 and 1 
orbital periods, and eccentricities between 0.999 and 0.9999, our 
improved guess is faster than converting to orbital elements and us¬ 
ing Danby’s by ^ 30%. In a manner analogous to that described in 
the previous section, we also solved Eq. [^perturbatively around 
X = dt/a = j3 dt/M in the regime » 1, but we found the 
second-order solution to be a slower initial guess than the simple 
X=p dt/M. 


3.5 Implementation of c-functions 

Einding a solution to Kepler’s equation is done iteratively and is 
thus the most expensive step in solving the Kepler problem. The 
iteration itself involves the calculation of multiple G-functions, 
which in turn require the calculation of c-functions. Thus, it is par¬ 
ticularly important to optimize these functions for both speed and 
accuracy. When calculating chaos indicators, we need cq, ci, C 2 , C 3 , 
C 4 and C 5 . If we are not integrating the variational equations, we 
only need cq, ci, C 2 and C 3 . 

We first ensure that z is smaller than 0.1 to guarantee that the 
series expansion of c in Eq. [^converges. We do this by dividing z 
repeatedly by 4. Note that divisions by powers of 2 are fast and ex¬ 
act in floating-point arithmetic. To calculate the series expansion, 
we need an inverse factorial for every term. Calculating this in¬ 
verse factorial by multiplying fioating-point numbers and then im¬ 
plementing a floating-point division would be very slow. We found 
that the fastest way to calculate the inverse factorial is to use a sim¬ 
ple lookup table. We checked that the series expansions of the c- 
functions converge very quickly for small z and thus we only store 
inverse factorials up to 1/34! in the lookup table. Any larger facto¬ 
rial would contribute less than one part in 10 ^^ to the sum and can 
thus be neglected (as we work in double fioating-point precision). 

We always calculate the first two terms in the series expan¬ 
sion. We then enter a loop and add more terms until the result no 
longer changes. Because z is small and the inverse factorials de¬ 
crease quickly, we are assured that the series will converge to a sin¬ 
gle floating-point number. This allows us to simply check whether 
the value changes from one iteration to the next, which is much 
faster than evaluating relative changes (cf. Sect. |3.2| ). 

Once the c-functions are calculated for the small z value, we 
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use the relations in Eqs. |20|21| with Eq. to calculate the c- 
functions for the original z value. 

Because this algorithm is an integral part of the integrator, we 
list the function to calculate c(z) in pseudo code: 


w <— 0 

while z > 0.1 do 

z <— z/4 
n n + I 

^ T\ ~ ^' T\ 
z <-z 

P^Z 

k^S 


repeat 


^ 4 ,prev ^ ^4 

p-z 

C4 C4 + p • ^ 

k k + \ 

C5 «- C5 + p • ^ 

k k + \ 


> Counter for quarter-angle formula. 

> Ensure that z is small. 


> Hard coded first two terms for C 4 . 

> p will the {-zy factor in the loop. 
> Third term in C 4 contains factor ^. 

> Keep old value to check for convergence. 

> 1 /^! comes from lookup table. 


until C4 — ^^4,prev 
C3 ^ I - Z • C5 

C2 <- I - Z • Q 
Cl <- 1 - z • C3 

while /I > 0 do 

z <— 4 • z 


> Converged? 
> Use Eq. 18 to get C 3 , C 2 and ci. 


> Apply quarter angle formula n times. 


C5 <— • (C5 -h C4 -h C3 -1- C2) 

Q «- I -cs •(! - Cl) 

C 3 g - z ■ C 5 
C2 «- 5 - z • C4 
Cl «- 1 - Z • C3 


n n - I 


Co ^ 1 - z • C 2 


3.7 A full integration in Jacobi coordinates 

The algorithms to convert to and from Jacobi coordinates that we 
describe in Sect. |3.1| are unbiased and fast. Nevertheless, we aim to 
avoid as many conversion as possible. 

As it turns out, we can reduce the number of conversions per 
timestep to two, one for the positions from Jacobi coordinates to the 
inertial frame, and one for the accelerations from the inertial frame 
to Jacobi accelerations. But note that this is only possible under 
the following assumptions: 1 ) the particle position and velocities 
are not changed in-between timesteps, e.g. manually by the user 
or by collisions, 2) outputs are not required at every timestep, 3) 
variational equations are not integrated, 4) no additional velocity- 
dependent forces are present. In such a case, an integration starting 
from an arbitrary inertial frame is achieved as follows: 

calculate Jacobi coordinates 

drift all particles under 7/Kepier for half a timestep, dt/2 
while t < ^max do 

calculate 1st part of //interaction in Jacobi coordinates 
update positions in the inertial frame 
calculate 2 nd part of //interaction in inertial frame 
convert accelerations from 2nd part to Jacobi accelerations 
apply kick from Jacobi accelerations to Jacobi velocities 
if not last timestep then 

drift all particles under //Kepler for a full timestep dt 

drift all particles under //Kepler for half a timestep dt/l 
update both positions and velocities in the inertial frame. 

Note that we never update the velocities in the inertial frame until 
the end of the simulation (or when an output is needed). We only 
convert the positions and velocities to Jacobi coordinates at the very 
beginning and not at every timestep. Besides the obvious speed¬ 
up, avoiding to go back and fourth between different coordinate 
systems reduces the build-up of round-off errors and thus makes 
the integrator more robust. 


3.6 Implementation of GauB / and g-functions 

The precise implementation of Gauss / and g functions matters for 
long term integrations. The straightforward implementation follow¬ 
ing ([Mikkola & Innan^|1999| ) leads to the / and g-functions in 
Eq. |27[ Note that for timesteps smaller than half an orbital period, 
the term MG 2 lro in / is small compared to the first term (which is 
just 1). The same argument holds true for g. We can define new / 


and g-functions 




/= 


MGi 


f = 


(37) 


TqT 


g = dt- MGi 


MG 2 

(38) 

1 = 



This allows us to rewrite the last step in solving the Kepler problem 
as 

r = [fro + «Vo) + ro v = (/ro + |vo) + Vq. (39) 

Although this step is algebraically equivalent to the original Eq.[^ 
we achieve higher precision. The reason is that we can now ensure 
that the small quantities in brackets are summed before they are 
added to the larger quantity (the initial value). We implement the 
same trick for the variational equations, Eqs.[^and[3T] 


3.8 LCN calculation 

To calculate the Lyapunov characteristic number and the Lyapunov 
timescale we need to perform a linear least square fit to the function 
Y(t). Thus we need the mean and the covariance of Y(t). Storing all 
previous values of Y(t) just to calculate its mean and covariance 
is inefficient. We therefore implement an efficient one-pass method 
described by |Pebay| < |2008| ). This method lets us calculate the LCN 
at every timestep in 0 (1) and has the further advantage of being 
numerically more robust than the standard implementation. 


4 NUMERICAL RESULTS 

In this section, we test the speed, accuracy and numerical stability 
of WHFast and compare it to other publicly available and widely 
used integrators. We begin by briefly defining our nomenclature 
for these other integrators and summarizing their properties. 

MERCURY is a mixed-variable symplectic integrator imple¬ 
mented in fortran and provided by the MERCURY package ( |Cham-| 
|bers & Migliorini||1997| ). This Wisdom-Holman style integrator 
uses high-order symplectic correctors. We directly call the fortran 
code without any modifications. 
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timestep, \oqlO{dt/t^^^) 


timestep, \oqlO{dt/t^^^) 


Figure 1. Tests of the Kepler-solver. Simulation with two bodies, integrated for 100 orbits with varying eccentricity and timestep. Left column: results using 
the standard WH integrator. Right column: results using our new WHFast integrator. Top row: relative energy error at the end of the simulation. Middle row: 
sign of the energy error at the end of the simulation. Bottom row: average runtime for one timestep. 
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SWIFTER-WHM is again a classical 2nd-order Wisdom-Holman 
integrator without symplectic correctors ( [Wisdom & Holni^ 
|1991| l. We use the integrator provided by the SWIFTER package. 
It is implemented in fortran and we directly call the SWIFTER exe¬ 
cutable without any modifications. 

SWIFTER-HELIO is also 2nd-order symplectic integrator with¬ 
out symplectic correctors ( [Duncan et al.|199^ . It uses democratic 
heliocentric coordinates. We again use the integrator provided by 
the SWIFTER package. It is implemented in fortran and we directly 
call the SWIFTER executable without any modifications. 

SWIFTER-TU4 is a 4th-order symplectic integrator. It is not 
a Wisdom-Holman integrator but splits the Hamiltonian in kinetic 
and potential terms ( [Gladman et al.[199T[ l. We also use the integra¬ 
tor provided by the SWIFTER package. It is implemented in fortran 
and we directly call the SWIFTER executable without any modifica¬ 
tions. 

For a more direct comparison, we also make use of an integra¬ 
tor that we simply refer to as WH. It is based on the SWIFTER-WHM 
integrator in SWIFTER but ported to C and available in the REBOUND 
( [Rein & Liu[2012[ l package. Like the SWIFTER-WHM integrator, it is 
a symplectic integrator that works in the heliocentric frame, and 
does not implement any symplectic correctors. Note that this is not 
the original integrator used by [Wisdom & Holman[ ( [199l] l, which is 
not publicly available. 

WHFast is C99 compliant. The C99 standard guarantees that 
fioating point operations are not re-ordered by the compiler (un¬ 
less one of the fast-math options is turned on). Because of that, 
the final positions and velocities of particles agree down to the last 
bit across different platforms. This makes WHFast platform inde¬ 
pendent and the simulation results reproducible. We verified this 
on different architectures (Linux, MacOSX), different CPUs (Intel 
Core i5-3427U, Intel Xeon E5-2697 v2, Intel Xeon E5-2620 v3) 
and different compilers (Apple LLVM 6.1.0, gcc 4.4.7). 

4.1 Two-body Kepler Solver 

The kernel of every Wisdom-Holman integrator is the Kepler 
solver. We describe our implementation in detail in Sections [3.2f 
[3.6[ Here, we test the Kepler solver using a two-body problem. 
The two body problem is invariant with respect to rescaling of 
the total mass, the mass ratio, the value of the gravitational con¬ 
stant and the orbital period. What does matter is the eccentric¬ 
ity of the orbit and the ratio of the timestep to the orbital period. 
We thus scan the parameter space in those two dimensions by in¬ 
tegrating two bodies for 100 orbital periods. We explore an ex¬ 
tremely wide parameter space. The eccentricities range from zero 
to 0.999 999 99 = 1-10"^. The range of timesteps goes from 0.1% 
of the orbital period all the way up to one orbital period. 

Fig .[^shows the performance of WHFast (right column) com¬ 
pared to WH (left column). The top row shows the absolute value of 
the relative energy error at the end of the simulation. The middle 
row shows the sign of the energy error. The bottom row shows the 
average runtime for a single timestep. The vertical lines visible in 
the top row correspond to times tep resonances ([Wisdom & Holm^ 
[1992[[T^ma & Wisdom| 199^ [Rauch & Holman|1999[l. 

One can see that WHFast is significantly more accurate than 
the standard WH integrator for the most important parts of param¬ 
eter space (eccentricities less than ~ 0.99). The relative energy is 
conserved better by two to three orders of magnitude. Most im¬ 
portantly, note that the energy error in the standard WH integrator 



• WH SWIFTER-HELIO •WHFAST-NOCOR 

• SWIFTER-WHM •MERCURY •WHFAST 

• SWIFTER-TU 4 


Figure 2. Relative energy error in simulations of the outer Solar System 
after 1000 Jupiter orbits as a function of the number of steps per orbit. 

is biased over large regions of the parameter space (there are large 
blue and red areas in the second row). On the other hand, WHFast 
has a random energy error throughout the parameter space. Having 
a biased energy error will lead to a long-term linear growth of the 
energy error (see below). 

In the entire parameter space explored, WHFast requires less 
time to complete a timestep than WH. The speed-up is typically be¬ 
tween 20% and 100%. Eor the integrations performed in this sec¬ 
tion, we convert to and from Jacobi coordinates at every timestep 
to provide a fair comparison. Thus, the speed-up and the energy- 
conservation properties of WHFast are in fact even better than 
shown here in any actual production run (see Sect [3.7[ ). 

4.2 Short Term Energy Conservation 

To compare the accuracy of the different integrators in a realistic 
test case, we run simulations of the outer Solar System for one 
thousand Jupiter orbits (12 000 years). We include the Sun and four 
massive bodies with approximate initial conditions corresponding 
to those of Jupiter, Saturn, Uranus and Neptune. In each simulation 
the initial conditions and masses are randomly perturbed by 0.1%. 
In Eig.|^ we plot the relative energy errors at the end of the simu¬ 
lation as a function of the number of timesteps imposed per Jupiter 
orbit. 

One can see that all the integrators except SWIFTER-TU4 are 
second-order schemes. Eor timesteps between 20% and 0.1% of the 
orbital period of Jupiter (50 to 1000 timesteps per orbit), their error 
decreases quadratically with decreasing timesteps. This is the error 
term Ebound introduced in Sect. [2.7[ 

However, decreasing the timestep also increases the number of 
fioating point operations. There will therefore be a timestep value 
at which the numerical round-off error dominates over the error as¬ 
sociated with the symplectic method itself Ebound- For that reason 
we find that for small timesteps, less than 0.1% of the shortest or¬ 
bital period, the errors of all integrators rise instead of decreasing 
further. Thus there is an optimum timestep Jtopt that yields the min¬ 
imum energy error. This optimum timestep depends on the length 
of the integration and will be larger for longer simulations. 

In Eig. one can see that the errors of WH, SWIFTER-WHM, 
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SWIFTER-HELIO and MERCURY rise very rapidly after reaching 
Jtopt. scaling as at least dr^ for the first decade. 

The optimum timestep for WHFast is roughly 0.1% of the 
shortest orbital period. However, WHFast’s error grows much more 
slowly with decreasing timestep than that of the other second-order 
integrators. In fact, the error is dominated by £rand and thus fol¬ 
lows as the number of timesteps A^steps increases as ~ dr^ if 
we keep the total integration time constant. Thus the behaviour of 
WHFast in Fig.|^for small timesteps can be seen as the first indica¬ 
tion that WHFast follows Brouwer’s law (see Sects. |2.7| and |4.4| ). 

The SWIFTER-TU4 integrator is the only other integrator we 
tested that seems to follow Brouwer’s law, but it performs poorly 
at large timesteps. This is expected, since unlike the other integra¬ 
tors, SWIFTER-TU4 does not assume a Keplerian splitting and must 
therefore take smaller timesteps to accurately reproduce the orbital 
motions. 

Integrators with symplectic correctors, MERCURY and WHFast, 
perform significantly better for long timesteps. Their energy con¬ 
servation is three orders of magnitude better (£^bound is three orders 
of magnitude smaller) compared to integrators without symplec¬ 
tic correctors. This is due to the mass ratio of Jupiter and the Sun 
being roughly 10“^. The order of the symplectic corrector is not 
very important for relatively high mass ratios such as these, i.e. a 
fifth-order symplectic corrector performs as well as an llth-order 
one. For much smaller mass ratios (when the mass ratio is less than 
the timestep ratio), higher-order symplectic correctors are advanta¬ 
geous. 

Note that Jtopt for almost all of the integrators is 10“^ orbital 
periods of Jupiter, i.e. 4 days. This is significant because Mercury’s 
orbital period is 88 days. Thus if we included Mercury in our sim¬ 
ulation, we would be very restricted in our timestep choice. We 
need more than 20 timesteps {dt ^ 4 days) to resolve Mercury’s 
orbit accurately. However, if we choose choose a timestep smaller 
than 4 days, we start to accumulate errors in the outer Solar Sys¬ 
tem. It is worth reiterating that the simulations shown in Fig. I^all 
ran for only 1000 orbits. If we ran a longer simulation with the 
same timestep, we would have more timesteps and thus accumu¬ 
late more round-off errors by the end of the simulation. One can 
therefore reach better energy conservation with a longer timestep. 
In other words, Jtopt is larger for longer integration times. 


4.3 Speed Comparison 

We run the same simulations as in Sect. |4.2| to compare the speed 
of the different integrators. Fig. shows the relative energy error 
as a function of runtime. The results show that no matter what the 
desired energy error is, WHFast is the fastest integrator. In the large 
timestep limit, the speed-up compared to MERCURY is roughly a fac¬ 
tor of 5. 

In the small timestep limit, dt < dt^^u we can only compare 
WHFast to SWIFTER-TU4, as all the other integrators’ errors are 
significantly larger (by 4 to 5 orders of magnitude) due to numer¬ 
ical roundoff errors (see below). SWIFTER-TU4 is as fast for small 
timesteps as WHFast but, as noted above, is unsuitable for large 
timesteps since it is not a Wisdom-Holman integrator. It is only 
shown here as a comparison. 
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Figure 3. Relative energy error in simulations of the outer Solar System 
after 1000 Jupiter orbits as a function of run time. 


4.4 Long Term Energy Conservation 

Let us finally address the most important benchmark, the long term 
energy conservation properties of WHFast compared to other inte¬ 
grators in a real world test case. In this section we only study the 
energy error, but other conserved properties like the angular mo¬ 
mentum behave the same way. In Fig. we show the time evolu¬ 
tion of the relative energy error in a simulation of the outer Solar 
System. As in Sect |4.2| we include the Sun and four massive bod¬ 
ies with approximate initial conditions corresponding to those of 
Jupiter, Saturn, Uranus and Neptune. The timestep for all simula¬ 
tions is 1.5 days. Note that this timestep is smaller than what one 
would typically choose for this kind of integration. However, with 
a 1.5 day timestep we reach machine precision for integrators that 
use symplectic correctors, allowing us to better quantify the long¬ 
term behaviour of WHFast. All the effects we discuss here are also 
present in simulations with longer timesteps, but they would man¬ 
ifest themselves in the relative energy error at a later time. We run 
four simulations for each integrator and randomly perturb the ini¬ 
tial conditions and masses by 0.1% in each simulation. We plot the 
individual simulations as thin lines, and the average error as a bold 
line. 

The lower relative energy bound is set by machine precision 
for all integrators, roughly 10"^^. WHFast and MERCURY, the inte¬ 
grators in our sample that have symplectic correctors, almost reach 
this limit early on in the simulation. The bound energy error Ebound 
is approximately 10"^"^. The integrators without symplectic correc¬ 
tors, SWIFTER-WHM and WH have an energy error roughly three order 
of magnitudes higher Ebound ~ 10"^^-^. 

From Fig. 1^ it is clear that the integrators MERCURY, WH and 
SWIFTER-WHM show a linear behaviour in the energy error at late 
times. This is due to the term Ebias- The Ebias term already domi¬ 
nates at early times (after 100 Jupiter orbits) for MERCURY because 
the symplectic correctors lower the value of Ebound- For WH and 
SWIFTER-WHM the Ebias term dominates after 10000 Jupiter orbits. 
This result shows that one or more steps in these integration algo¬ 
rithms are biased. We found that the two main contributions were 
the inaccurate implementation of the rootfinder for Kepler’s equa¬ 
tion and the conversions to and from Jacobi coordinates. In WHFast, 
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-MERCURY - WH -SWIFTER-WHM 

- WHFAST-NOCOR - WHFAST 


Figure 4. Relative energy error in simulations of the outer Solar System as a function of time for different symplectic integrators. The timestep for all 
simulations is 1.5 days. 


£bias is absent, showing that its implementation is completely unbi¬ 
ased. 

Since all integrators are implemented in double floating-point 
precision and use the same timestep, they all have roughly the same 
error term £^rand- However, it is only visible in Fig.j^for the WHFast 
integrator. For all other integrators the linearly growing term £^bias 
dominates over £^rand- 

If we increase the timestep, the linear error growth will show 
up at a later time because £^bound will be larger. However, it is still 
present at all times. Let us think of a symplectic integrator as an ex¬ 
act integrator for a perturbed Hamiltonian Ti with high frequency 
terms added compared to in Eq.[^ see e.g. [Wisdom et al.| { |1996] ). 
Then the quantity related to the energy error for TY, let us call this 
E, should be conserved exactly at all times (that is the idea of a 
symplectic integrator). However, if the implementation is biased, 
E will undergo a linear growth at all times. With WHFast, we im¬ 
prove the conservation of E by many orders of magnitude in any 
integration, regardless of timestep. 

This difference could have important implication for the dy¬ 
namical evolution of the system and could for example push it from 
a stable to an unstable region of parameter space. We plan to study 
the effect of different integrators on systems near a chaotic/non- 
chaotic separatrix in a follow up paper. 


5 CONCLUSIONS 

In this paper, we presented WHFast, a new implementation of a 
symplectic Wisdom-Holman integrator. Key advantages and im¬ 
provements over other publicly available implementations of sym¬ 
plectic integrators are: 

WHFast is faster by a factor of 1.5 to 5. Of that, a 50% speedup 
comes from the improved Kepler solver, where we use a fast con¬ 
vergence criteria for Newton’s method and an efficient implemen¬ 
tation of c and G-functions. The remainder of the speedup is due to 
combining drift steps at the end and beginning of each timestep and 
to only converting to and from Jacobi coordinates when needed. 


The Kepler solver is more accurate and unbiased. We achieve 
this thanks to improvements to the convergence criteria in Newton’s 
method, a Laguerre-Conway solver for highly eccentric orbits with 
long timesteps, the high accuracy implementations of the c and G- 
functions and a careful ordering of floating-point operations. 

We remove the secular energy error that grows linearly with 
integration time. This is due to two improvements. First, the un¬ 
biased Kepler solver. Second, the improved and also unbiased co¬ 
ordinate transformations to and from Jacobi coordinates. To our 
knowledge, WHFast is the first publicly available implementation 
of a Wisdom-Holman integrator that follows Brouwer’s law over 
long timescales for small enough timesteps and does not show a 
linear growth in the energy error. 

We implement variational equations that allow us to compute 
the Lyapunov timescale and the MEGNO. Our algorithm to calcu¬ 
late the Lyapunov timescale uses a numerically stable algorithm 
that is based on a one-pass covariance Alter. The variational equa¬ 
tions do not require us to solve Kepler’s equation and are thus very 
inexpensive to calculate. 

Symplectic correctors of order 3, 5, 7, and 11 are imple¬ 
mented. These symplectic corrector allow for high-accuracy simu¬ 
lations of systems with small mass ratios. Even for relatively mas¬ 
sive planets like those in the Solar System, symplectic correctors 
achieve an improvement of three orders of magnitude. Eor long in¬ 
tegrations, the performance cost of symplectic correctors is negli¬ 
gible and so our default setting uses an llth-order corrector. 

WHFast lets the centre-of-mass move freely during an integra¬ 
tion. We integrate an additional degree of freedom in order for our 
integrator to work in any inertial frame, i.e. one is not restricted to 
the heliocentric or barycentric frame. Additionally, we do not tie 
our implementation to a specific choice of units. 

The integrator is available as an easy to use python module. 
The module works on both python 2 and 3. It can be installed on 
most Unix and MacOS systems with a single command: 

pip install rebound 

The following python script imports the rebound module, adds par- 
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tides to the simulation, selects an integrator and timestep and runs 
the integration. 

import rebound 

rebound.add(m=l) 

rebound.add(m=0.001, a=l.) 

rebound.add(m=0.001, a=2., e=0.1) 

rebound.integrator = ’whfast’ 

rebound.dt = 0.01 

rebound.integrate(6.2831) 

More complicated examples and the source code of WHFast (writ¬ 
ten in C, compliant with the C99 standard) can be found in the 
REBOUND package. REBOUND includes several other integrators, col¬ 
lision detection algorithms, a gravity tree code and much more. 
The REBOUND git repository is hosted at https://github.com/ 
hannorein/rebound 

We also provide an experimental hybrid integrator for sim¬ 
ulations in which close encounters occur. The hybrid integrator 
switches over to a high-order non-symplectic integrator (IAS 15, 
|Rein & Spiegel|20T5] ) during a close encounter. A detailed discus¬ 
sion of this integrator and its properties will be given in a follow-up 
paper. 

We hope that with the speed and accuracy improvements, 
WHFast will become the go-to integrator package for short and 
long-term orbit simulations of planetary systems. 
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